# -*- coding: utf-8 -*-
__author__ = 'Dandelion'
__email__ = 'xushiluo@163.com'

import os
import subprocess
import ConfigParser
import csv
from qgis.core import *
from qgis.utils import iface
from osgeo import gdal


class LandClassifier:

    #从ini文件读取OTB路径的设置
    @staticmethod
    def readIniItem(itemKey, iniFileName="otbSettings.ini"):
        if not os.path.exists(iniFileName):
            return None
        config = ConfigParser.ConfigParser()
        config.readfp(open(iniFileName))
        return config.get("Default", itemKey)

    @staticmethod
    def findOtbPath():
        folder = None
        testfolder = os.path.join(os.path.dirname(QgsApplication.prefixPath()),
                                  os.pardir, "bin")
        if os.path.exists(os.path.join(testfolder, "otbcli.bat")):
            folder = testfolder
        #如果folder是None，那么从ini文件读取otbFolder的路径
        if folder is None:
            folder = LandClassifier.readIniItem("otbPath")
        #如果从ini文件读取也未找到otbPath，那么将其设置为c:\OSGeo4W\bin
        if folder is None:
            folder = "c:\\OSGeo4W\\bin"
        return folder

    @staticmethod
    def findOtbLibPath():
        folder = None
        testfolder = os.path.join(os.path.dirname(QgsApplication.prefixPath()), "orfeotoolbox", "applications")
        if os.path.exists(testfolder):
            folder = testfolder
        #如果folder是None，那么从ini文件读取otbFolder的路径
        if folder is None:
            folder = LandClassifier.readIniItem("otbLibPath")
        if folder is None:
            folder = "c:\\OSGeo4W\\apps\\orfeotoolbox\\applications"
        return folder

    #计算图像的全局均值和各波段的标准差
    @staticmethod
    def execComputeImagesStatistics(inputImageFullPath, outputName=None):
        if outputName is None:
            return None
        os.putenv('ITK_AUTOLOAD_PATH', LandClassifier.findOtbLibPath() )
        cliName = "otbcli_ComputeImagesStatistics"
        path = LandClassifier.findOtbPath()
        commands = []
        commands.append(path + os.sep + cliName)
        commands.append("-il")
        commands.append(inputImageFullPath)
        commands.append("-bv")
        commands.append(0)
        commands.append("-out")
        commands.append(outputName)
        fused_command = ''.join(['"%s" ' % c for c in commands])
        proc = subprocess.Popen(fused_command, shell=True, stdout=subprocess.PIPE, stdin=subprocess.PIPE,stderr=subprocess.STDOUT, universal_newlines=True).stdout
        loglines = []
        loglines.append("ComputeImagesStatistics output")
        for line in iter(proc.readline, ""):
            loglines.append(line)
        return "".join(loglines)

    # 构建Bayes分类器
    @staticmethod
    def execBayesClassfier( inputImageFullPath, inputVectorFileName, statisticsXmlName, bayesModelName=None, bayesConfusionMatrixName=None, classFieldName="class"):
        """
        构建一个贝叶斯分类器，用于后续分类
        :param inputImageFullPath:输入的图像
        :param inputVectorFileName: 感兴趣区域（ROI）shape文件
        :param statisticsXmlName: 统计信息xml文件
        :param bayesModelName: 输出的分类器完整文件名，后缀txt
        :param bayesConfusionMatrixName: 输出的分类器混淆矩阵文件名，后缀csv
        :return:字符串，执行信息
        """
        os.putenv('ITK_AUTOLOAD_PATH', LandClassifier.findOtbLibPath() )
        cliName = "otbcli_TrainImagesClassifier"
        path = LandClassifier.findOtbPath()
        commands = []
        commands.append(path + os.sep + cliName)

        commands.append("-io.il")
        commands.append(inputImageFullPath)
        commands.append("-io.vd")
        commands.append(inputVectorFileName)
        commands.append("-io.imstat")
        commands.append(statisticsXmlName)

        commands.append("-elev.default")
        commands.append(0)
        commands.append("-sample.mt")
        commands.append(1000)
        commands.append("-sample.mv")
        commands.append(1000)
        commands.append("-sample.vtr")
        commands.append(0.5)
        commands.append("-sample.vfn")
        commands.append(classFieldName)
        commands.append("-classifier")
        commands.append("bayes")
        commands.append("-rand")
        commands.append(0)

        commands.append("-io.confmatout")
        commands.append(bayesConfusionMatrixName)

        commands.append("-io.out")
        commands.append(bayesModelName)

        fused_command = ''.join(['"%s" ' % c for c in commands])
        proc = subprocess.Popen(fused_command, shell=True, stdout=subprocess.PIPE, stdin=subprocess.PIPE,stderr=subprocess.STDOUT, universal_newlines=True).stdout
        loglines = []
        loglines.append("build Bayes Classfier output:")
        for line in iter(proc.readline, ""):
            loglines.append(line)
        return "".join(loglines)

    #执行分类操作
    @staticmethod
    def ExecClassify( inputImageFullPath,statisticsXmlName,bayesModelName, reusltName=None ):
        """
        根据生成的分类器执行分类操作
        :param inputImageFullPath:待分类的图像
        :param statisticsXmlName: 待分类图像的统计信息
        :param bayesModelName: 分类器的完整文件名，是一个txt文件
        :param reusltName: 分类结果图像的完整文件名
        :return:执行情况字符串
        """
        if reusltName is None:
            return None
        cliName = "otbcli_ImageClassifier"
        path = LandClassifier.findOtbPath()
        commands = []
        commands.append(path + os.sep + cliName)

        commands.append("-in")
        commands.append(inputImageFullPath)
        commands.append("-model")
        commands.append(bayesModelName)
        commands.append("-imstat")
        commands.append(statisticsXmlName)
        commands.append("-ram")
        commands.append(128)

        commands.append("-out")
        commands.append(reusltName)

        fused_command = ''.join(['"%s" ' % c for c in commands])
        proc = subprocess.Popen(fused_command, shell=True, stdout=subprocess.PIPE, stdin=subprocess.PIPE,stderr=subprocess.STDOUT, universal_newlines=True).stdout
        loglines = []
        loglines.append("Classification output:")
        for line in iter(proc.readline, ""):
            loglines.append(line)
        return "".join(loglines)

    #执行颜色映射操作
    @staticmethod
    def ExecColorMapping( classifyImageFullPath, reusltName=None ):
        if reusltName is None:
            return None
        os.putenv('ITK_AUTOLOAD_PATH', LandClassifier.findOtbLibPath())
        cliName = "otbcli_ColorMapping"
        path = LandClassifier.findOtbPath()
        commands = []
        commands.append(path + os.sep + cliName)

        commands.append("-in")
        commands.append(classifyImageFullPath)
        commands.append("-ram")
        commands.append("128")
        commands.append("-op")
        commands.append("labeltocolor")
        commands.append("-method")
        commands.append("continuous")
        commands.append("-method.continuous.lut")
        commands.append("relief")
        commands.append("-method.continuous.min")
        commands.append("1")
        commands.append("-method.continuous.max")
        commands.append("5")

        commands.append("-out")
        commands.append(reusltName + u" uint8")
        fused_command = ''.join(['"%s" ' % c for c in commands])
        proc = subprocess.Popen(fused_command, shell=True, stdout=subprocess.PIPE, stdin=subprocess.PIPE,stderr=subprocess.STDOUT, universal_newlines=True).stdout
        loglines = []
        loglines.append("colormapping:")
        for line in iter(proc.readline, ""):
            loglines.append(line)
        return "".join(loglines)

    # 波段运算
    @staticmethod
    def CalcBandMath(inputFileNameList, expression, resultName):
        os.putenv('ITK_AUTOLOAD_PATH', LandClassifier.findOtbLibPath())
        if resultName is None:
            return None
        cliName = "otbcli_BandMath"
        path = LandClassifier.findOtbPath()
        commands = []
        commands.append(path + os.sep + cliName)

        inputFullPathFileName = " ".join(inputFileNameList)
        commands.append("-il")
        commands.append(inputFullPathFileName)
        commands.append("-ram")
        commands.append("128")
        commands.append("-exp")
        commands.append(expression)

        commands.append("-out")
        commands.append(resultName)

        fused_command = ''.join(['"%s" ' % c for c in commands])
        proc = subprocess.Popen(fused_command, shell=True, stdout=subprocess.PIPE, stdin=subprocess.PIPE,stderr=subprocess.STDOUT, universal_newlines=True).stdout
        loglines = []
        loglines.append("Band Math for change detection:")
        for line in iter(proc.readline, ""):
            loglines.append(line)
        return ''.join(loglines)

    # -------------------------------------------------------------------------------------------------
    @staticmethod
    def getComputeImagesStatisticsFileName(inFileName):
        """ 获取输入文件统计信息的完整文件名 """
        fpathandname, fext = os.path.splitext(inFileName)
        outputName = fpathandname + "_ComputeImagesStatistics.xml"
        return outputName

    @staticmethod
    def getClassifierModelName(inFileName,classifierName="bayesClassifier"):
        """
            获取输入文件的分类器模型完整文件名
        """
        fpathandname, fext = os.path.splitext(inFileName)
        outputName = fpathandname + "_" + classifierName + ".txt"
        return  outputName

    @staticmethod
    def getConfusionMatrixName(inFileName, conMatName="bayesConfusionMatrix"):
        """
        获取输入文件的混淆矩阵完整文件名
        """
        fpathandname, fext = os.path.splitext(inFileName)
        outputName = fpathandname + "_" + conMatName + ".csv"
        return outputName

    @staticmethod
    def getClassifyResultName(inFileName, classifierName="bayesResult"):
        """
        获取输入文件分类后的完整文件名
        """
        fpathandname, fext = os.path.splitext(inFileName)
        outputName = fpathandname + "_" + classifierName + fext
        return outputName

    @staticmethod
    def getColorMappingName(inFileName, colorTable="relief"):
        """
        获取分类后结果作颜色映射后的完整文件名
        """
        fpathandname, fext = os.path.splitext(inFileName)
        outputName = fpathandname + "_colorMapping_" + colorTable + fext
        return  outputName

    @staticmethod
    def getChangeDetectionName(baseInputFileName,secondInputFileName):
        path,baseImageFileName=os.path.split(baseInputFileName)
        baseImageBaseName, fext = os.path.splitext(baseImageFileName)
        secondImageFileName = os.path.basename(secondInputFileName)
        secondImageBaseName, fext2 = os.path.splitext(secondImageFileName)
        outputName = path + "/" + baseImageBaseName + "_" + secondImageBaseName + "_changeDetection" + fext
        return outputName

    # -------------------------------------------------------------------------------------------------

    #获取所有的栅格图层
    @staticmethod
    def getRasterLayers():
        layers = QgsMapLayerRegistry.instance().mapLayers().values()
        raster = []
        for layer in layers:
            if layer.type() == layer.RasterLayer:
                if layer.providerType() == 'gdal':  # only gdal file-based layers
                    raster.append(layer)
        return raster

    #获取所有矢量图层
    @staticmethod
    def getVectorLayers(shapetype=[-1]):
        layers = QgsMapLayerRegistry.instance().mapLayers().values()
        vector = []
        for layer in layers:
            if layer.type() == layer.VectorLayer:
                if shapetype == [-1] or layer.geometryType() in shapetype:
                    uri = unicode(layer.source())
                    if not uri.lower().endswith('csv') \
                            and not uri.lower().endswith('dbf'):
                        vector.append(layer)
        return vector

    #获取图层扩展的文件名
    @staticmethod
    def getExtendedLayerName(layer):
        authid = layer.crs().authid()
        if authid is not None:
            return layer.name() + ' [' + str(authid) + ']'
        else:
            return layer.name()

    # 根据图层对象获取输入图层的路径
    @staticmethod
    def getInputLayerPath(layerObj):
        lyrPath = None
        if isinstance(layerObj, QgsRasterLayer):
            lyrPath = unicode(layerObj.dataProvider().dataSourceUri())
        elif isinstance(layerObj, QgsVectorLayer):
            lyrPath = unicode(layerObj.source())
        else:
            lyrPath = unicode(layerObj)
            layers = LandClassifier.getRasterLayers()
            for layer in layers:
                if layer.name() == lyrPath:
                    lyrPath = unicode(layer.dataProvider().dataSourceUri())
        return lyrPath

     # -------------------------------------------------------------------------------------------------
    #计算每个灰度值在图像中所占的百分比
    @staticmethod
    def execComputeStatistics(inputRasterFileName,outputCSVFileName, readBandNo=1):
        # Open the input raster file
        # register the gdal drivers
        gdal.AllRegister()

        # Open and assign the contents of the raster file to a dataset
        dataset = gdal.Open(inputRasterFileName, GA_ReadOnly)
        if dataset is None:
            return False
        #pydevd.settrace('localhost',  port=53100,  stdoutToServer=True,  stderrToServer=True)
        #只读取第readBandNo个波段
        oldArray = dataset.GetRasterBand(readBandNo).ReadAsArray(0, 0, dataset.RasterXSize, dataset.RasterYSize).astype(numpy.float64)

        sizeTuple = oldArray.shape
        if len(sizeTuple) == 3:
            bandNums=sizeTuple[0]   #number of bands
            imgHeight = sizeTuple[1]
            imgWidth = sizeTuple[2]
        elif len(sizeTuple) == 2:
            imgHeight = sizeTuple[0]
            imgWidth = sizeTuple[1]

        TotalPixelNum = imgHeight*imgWidth
        valueLabelDict={}   #保存值标签，例如一幅图像可能只有2, 4，8这三个灰度值
        for val in oldArray.flat:
            tempValue = int(val)
            if tempValue not in valueLabelDict:
                valueLabelDict[tempValue] = 0   #每个灰度值的像素点数初始为0

        # 统计每个每个灰度值的像素点数，例如灰度值为2的像素共有358个，那么字典valueLabelDict[2]的值为358
        for val in oldArray.flat:
            tempValue = int(val)
            valueLabelDict[tempValue] = valueLabelDict[tempValue] +1

        #统计各个灰度值的百分比
        percentLabelDict={}
        for k in valueLabelDict:
            percentLabelDict[k] = float(valueLabelDict[k]) * 100.0 / float(TotalPixelNum)

        csvfile = file(outputCSVFileName, 'wb')
        writer = csv.writer(csvfile)
        writer.writerow(list(valueLabelDict))

        data = [
            valueLabelDict.values(),
            percentLabelDict.values()
        ]
        writer.writerows(data)

        csvfile.close()